home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 November / macformat-018.iso / Utility Spectacular / Utilities / Calc / lib / bernoulli.cal next >
Encoding:
Text File  |  1992-02-24  |  1.4 KB  |  63 lines  |  [TEXT/????]
  1. /*
  2.  * Copyright (c) 1992 David I. Bell
  3.  * Calculate the Nth Bernoulli number B(n).
  4.  * This uses the following symbolic formula to calculate B(n):
  5.  *
  6.  *    (b+1)^(n+1) - b^(n+1) = 0
  7.  *
  8.  * where b is a dummy value, and each power b^i gets replaced by B(i).
  9.  * For example, for n = 3:
  10.  *    (b+1)^4 - b^4 = 0
  11.  *    b^4 + 4*b^3 + 6*b^2 + 4*b + 1 - b^4 = 0
  12.  *    4*b^3 + 6*b^2 + 4*b + 1 = 0
  13.  *    4*B(3) + 6*B(2) + 4*B(1) + 1 = 0
  14.  *    B(3) = -(6*B(2) + 4*B(1) + 1) / 4
  15.  *
  16.  * The combinatorial factors in the expansion of the above formula are
  17.  * calculated interatively, and we use the fact that B(2i+1) = 0 if i > 0.
  18.  * Since all previous B(n)'s are needed to calculate a particular B(n), all
  19.  * values obtained are saved in an array for ease in repeated calculations.
  20.  */
  21. global Bn, Bnmax;
  22. mat Bn[1001];
  23. Bnmax = 0;
  24.  
  25.  
  26. define B(n)
  27. {
  28.     local    nn, np1, i, sum, mulval, divval, combval;
  29.  
  30.     if (!isint(n) || (n < 0))
  31.         quit "Non-negative integer required for Bernoulli";
  32.  
  33.     if (n == 0)
  34.         return 1;
  35.     if (n == 1)
  36.         return -1/2;
  37.     if (isodd(n))
  38.         return 0;
  39.     if (n > 1000)
  40.         quit "Very large Bernoulli";
  41.  
  42.     if (n <= Bnmax)
  43.         return Bn[n];
  44.  
  45.     for (nn = Bnmax + 2; nn <= n; nn+=2) {
  46.         np1 = nn + 1;
  47.         mulval = np1;
  48.         divval = 1;
  49.         combval = 1;
  50.         sum = 1 - np1 / 2;
  51.         for (i = 2; i < np1; i+=2) {
  52.             combval = combval * mulval-- / divval++;
  53.             combval = combval * mulval-- / divval++;
  54.             sum += combval * Bn[i];
  55.         }
  56.         Bn[nn] = -sum / np1;
  57.     }
  58.     Bnmax = n;
  59.     return Bn[n];
  60. }
  61.  
  62. print 'B(n) defined';
  63.